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Abstract 

An improved numerical solver for the unified solution of compressible and incom- 
pressible fluids involving interfaces is proposed. The present method is based on 
the CIP-CUP (Cubic Interpolated Propagation / Combined, Unified Procedure) 
method, which is a pressure-based semi-implicit solver for the Euler equations of 
fluid flows. In Part I of this series of articles [M. Ida, Comput. Phys. Commun. 
132 (2000) 44], we proposed an improved scheme for the convection terms in the 
equations, which allowed us discontinuous descriptions of the density interface by 
replacing the cubic interpolation function used in the CIP scheme with a quadratic 
extrapolation function only around the interface. In this paper, as Part II of this 
series, the multi-time-step integration technique is adapted to the CIP-CUP integra- 
tion. Because the CIP-CUP treats different-nature components in the fluid equations 
separately, the adaptation of the technique is straightforward. This modification al- 
lows us flexible determinations of the time interval, which results in an efficient 
and accurate integration. Furthermore, some additional discussion on our methods 
is presented. Finally, the application results to composite flow problems such as 
compressible and incompressible Kelvin-Helmholtz instabilities and the dynamics 
of two acoustically coupled deformable bubbles in a viscous liquid are provided. 
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1 Introduction 



This series of articles presents an improved solver for a challenging problem, 
the unified solution of compressible and incompressible fluids. After Harlow 
and Amsden proposed the ICE algorithm as a fully implicit solver for fluid 

equations in a conservative form [1], some approaches for this purpose have 
been investigated, such as the pressure-based semi-implicit algorithms [2,3,4], 
approaches based on the asymptotic expansion with respect to the local Mach 
number [5,6], and the boundary condition capturing method that treats com- 
pressible and incompressible materials separately [7]. (See also recent reviews 
[8,9] for more details.) Among them, the CIP-CUP method, a pressure-based 
semi-implicit algorithm proposed by Yabe and Wang [3] , has already been ap- 
plied to many practical multi-material problems such as the laser machining of 
a metal plate [10], comet Shoemaker- Levy 9's colhsion with the planet Jupiter 
[11], the interaction of a shock wave and a liquid drop [12], and the milk-crown 
formation on a liquid surface [13], and has been proven to be an efficient and 
robust solver for the unified solution [14,15,16,17]. In this method, the Euler 
(or the Navier-Stokes) equations for fluid flows in a non-conservative form are 
selected as the governing equations, and the convection terms in the equa- 
tions are solved explicitly by the CIP method [18], while the acoustic terms 
are solved implicitly by the CUP method [3]. In Part I of this series [19], an 
improved solver for the convection terms was constructed using both an in- 
terpolation and an extrapolation function to describe the spatial profile of the 
density of materials. In the improved method, the cubic interpolation function 
used in the CIP scheme is replaced with the quadratic extrapolation functions 
constructed under some constraints for guaranteeing stability, only in a cell 
containing an interface between materials. This method allows us to solve the 
density interface with no dissipation across the interface, and is applicable to 
compressible flow problems, i.e., to variable-density problems. 

As Part II of this series, this paper presents some further improvements for 
the CIP-CUP method, and also gives some application results. In Sec. 3, we 
attempt to incorporate the concept of the multi-time-step (MTS) integration 
techniques into the CIP-CUP algorithm. The MTS integration techniques, 
which allow us the efficient integration of dynamical systems involving some 
different time scales, have been used to solve a variety of problems, such as 
gravitational A^-body problems [20,21], molecular dynamics [22,23,24], and at- 
mospheric fiow problems [25] . This technique reduces the computational effort 
by integrating only rapidly varying components of a system (e.g., strongly ac- 
celerated particles or waves having fast phase speeds) with a small time inter- 
val, and others with a larger time interval. In the present study, based on this 
concept, different-nature terms in the fluid equations (convection, acoustic, 
and other terms) are solved using different time intervals. This modification 
allows us a more flexible determination of the time interval, and improves the 



2 



accuracy and efficiency of tlic CIP-CUP. As will be shown, the adaptation 
of the MTS concept to the CIP-CUP integration is straightforward, because 
the CIP-CUP separately treats the terms having different natures, i.e., those 
having different characteristic time scales. 

Section 4 presents a modified scheme to update the spatial derivatives during 
the computational steps except for the convection parts. As was shown in 
Part I, the spatial derivatives of the dependent variables are used explicitly as 
additional dependent variables to solve the convection terms. In the convection 
process, the derivatives are updated to obey the equations derived by taking 
derivatives of the governing equations [18,19]. Those spatial derivatives should 
be updated by some way also in the non-convection processes. In Refs. [18,26], 
simple methods for this purpose have been proposed using classical centered 
differences. As was discussed in Part I, however, the use of such a classical 
discretization gives rise to the numerical dispersion and dissipation around 
the interfaces at which the spatial derivatives of the dependent variables tend 
to be discontinuous. We modify the centered scheme by adopting a simple 
extrapolation. 

In Sec. 5, the applicability of the present methods to the compressible Navier- 
Stokcs equations with a surface tension term is investigated, and in the Ap- 
pendixes, some additional discussions regarding the averaging of the density 
(needed to solve the terms other than the convection) and the treatment of 
the surface-tension term are given. Section 6 presents some application results 
for the compressible and incompressible Kelvin-Helmholtz instabilities and 
the multibubble dynamics in an acoustic field, and Sec. 7 presents concluding 
remarks. 



2 The CUP method and its variants 

2.1 The CUP method 

The CIP-CUP method [3] is a semi-imphcit solver for the Euler equations of 
fluid flows: 



-pV • u. 





Vp F 



P P 



(2) 




-pC|V • u. 



(3) 
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where p, u, p, and Cs denote the density, the velocity vector, the pressure, 
and the local sound speed, respectively, and F may contain the viscosity, the 
surface tension, and external forces. This method separately solves the terms in 
Eqs. (l)-(3) of different natures by a time splitting technique. The convection 
parts of these equations, 

| + u.V, = 0. 

du „ „ 
^ + u.Vu = 0. 

dp „ 



are solved by the CIP method [18] (In the present work, of course, the hybrid 
interpolation-extrapolation method proposed in Part I is adapted), while the 
acoustic parts, 

|^-.V.u. (4) 



du Vp 

'dt " ~V 



(5) 



I = -pC|V • u, (6) 



are solved by the CUP method, which is an implicit finite difference method. 
The remaining part, 

- = - (7) 



may be solved by some existing methods, such as the finite difference or finite 
volume method. (We refer to Eq. (7) as the "additional part".) Solving these 
parts successively completes one step of the CIP-CUP time integration. In the 
following, the CUP method for the acoustic parts is concretely reviewed. 

Discretizing the time derivatives on the LHS of Eqs. (4)-(6) and estimating 
the spatial derivatives on the RHS with future values, one obtains 

^ ^ ^-p-V-u"+\ (8) 
At ~ ^ ^ 
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- -P*Q^V • u"+S (10) 



where n is the number of the time step, At is the time interval, and the 
quantities with the superscript * indicate the values after solving the parts 
other than the acoustic parts. (The acoustic parts should be solved at the final 
stage of a time step, because, as is shown below, solving these parts enforces 
the divergence-free condition for an incompressible fluid.) Taking divergence of 
Eq. (9) and substituting it into Eq. (10) yield the following pressure equation: 

„n+l _ „* ^„n+l 

^ , ^ = p*Cf/\t V ■ — p*CfV ■ u*. (11) 



(Substituting Eq. (6) into this and rewriting as 



finds that Eq. (11) corresponds to a first-order approximation of the wave 
equation in terms of the pressure.) Almost the same pressure equation is given 
in Ref. [4], but discretization is done by a finite element technique. 

After solving Eq. (11) to get p""^"^, the velocity is updated exphcitly by Eq. 
(9). Also, the density is updated exphcitly by 



= + ^ . (12) 



given by Eqs. (8) and (10). Generally, the spatial discretization for the above 
equations is performed using the 2nd-order centered finite differencing on the 
staggered grids. 

When the sound speed is infinite, the pressure equation (11) is reduced to an 
elliptic equation, 

V.^ = ^. (13) 

p* At ^ ^ 



which corresponds to that used in the SMAC algorithm [27] for incompress- 
ible fiows. This result reveals that the CUP method is applicable to both 
compressible and incompressible fiuids. 

In the case where the CIP or its variant is used to solve the convection parts, 
we also need to update the spatial derivatives, used explicitly in those schemes. 
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in the non-convection processes. Simple schemes for this purpose have already 
been proposed [18,26]. For a two-dimensional case, the scheme is represented 
as 



fn+l xn+1 X* £* 



fn+1 fn+1 f* £* 

Q rn+l f) f* _ Ji,3+'^ -/hj-^ Ji,j-1 p-\ 



where / indicates an dependent variable, d^f = df/dx, dyf = df/dy, the 
subscripts i and j indicate {x, y) = {i h, j h), and h is the grid spacing (assumed 
to be constant for simplicity). In the case where the cross derivative d^yf is 
used in the convection process [26,19], the following equation is additionally 
used: 



£n+l fn+1 fn+1 , fn+1 

r) f"+l _ f * _ Ji+lJ+l Ji-l,j+l Ji+l,j-l ~^ Ji-l,j-l 

f* f* f* I f* 

Ji+l,j+l Ji-l,j+l Ji+l,j-l "I" Ji-l,j-l 



(16) 



Equations (14)-(16) can be solved explicitly using the quantities obtained 
before and after solving the non-convection parts. 



2.2 Improved variants of the pressure equation 



In 1994 [28], Ito proposed an improved variant of the pressure equation [Eq. 
(11)] by incorporating the concept of the exponential method for a heat- 
conduction equation [29]. The 1-D formula of the variant is represented by 

= P*Cf^t [a [P-y^ -f (1 - a) ] - (17) 



where a is the weighting factor determined theoretically as 

a{E) ^ ^ 



l-exp(-£;) E' 
, 1 1 \ * fCs*At^ ^ 



Pi+l/2 Pi-1/2 



h 



and Pi±i/2 is the density sX x = h{i + 1/2), determined approximately with p* 
and p*^i (see Appendix A). For — > oo (i.e., Csl^t/h — > cxd), a converges to 
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1 and Eq. (17) is reduced to the elliptic equation (13); therefore, this variant 
is applicable to incompressible flows. For E ^ 0, a becomes 1/2, resulting in 
a higher resolution than that of the conventional CUP [28,30]. We adopted 
this variant in the present study. 

The 2-D formula of this variant (not shown in Ito's paper) may be represented 

by 



„n+l _ „* _ _ 



At 



where 



Gl = a{El)(P^\ +(l-a(El))(^j , 



G2^a{E2)(^^^ + (1 - a(£;2)) (^pj , 



1_ ]_\ . fCslAt' 

Pi+l/2,j Pi-l/2,j 



2 



Ei={- — + - — \PlA^ 



E2 = + 



1 1 \ , fCsl^At^' 



f 1,3+1/2 Pi,j-i/2) \ J 
We simplify this as follows to reduce the computational efforts: 



P"^' P* ^p*CfAtHE,^,^)V-^^ + {l-a{E^,^))V 



At p* p* 

-p*C*s^V ■ u*, (18) 

where 

£;^ax = max(El,E2). 

This formula requires only one weighting factor. This simplification may be 
valid because El a; E2 generally. 

In Ref. [30], an alternative variant was proposed. Though this variant can 
provide higher resolution of the sound wave than that obtained by Ito's variant, 
it is likely that the use of this variant as a part of the solver for fluid flows 
sometimes causes great violation of mass conservation [31]. We therefore did 
not adapt it. 
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3 Multi-time-step integration 



In many fluid flow problems, the time scale of a signiflcant phenomenon is de- 
termined by the flow velocity (i.e., the characteristic speed of the convection 
parts), rather than by the sound speed or other characteristic speed. In certain 
cases, however, the acoustic parts (or others) have much smaller time scales 
(i.e., much greater characteristic speeds), resulting in an excessive restriction 
on the time interval. At. The CIP-CUP overcomes this problem by solving 
the acoustic parts implicitly. Such an approach, however, guarantees only the 
stability, and probably provides inaccurate results in certain situations. For 
example, under an intermediate condition of flows such as weakly compress- 
ible flows, the acoustic parts may need to be solved by using a time interval 
of sufficiently small kc that makes Ku very small, since the compressibility 
described by the acoustic parts plays an important role, where 

Kc = max(Cs')At//i and Ku = max(|u|)At//i. 



To achieve an efficient computation even in such a situation, we adopt the 
multi-time-step (MTS) integration technique to the CIP-CUP method. 

The adaptation of the MTS integration to the CIP-CUP is straightforward 
because the components of different time scales (convection, sound propaga- 
tion, and others) are treated separately. The time propagation performed in 
the CIP-CUP can be represented schematically as follows: 

U** = Ll(At)U", (19) 

U* = L2{At) U**, (20) 

U"+^ = L3{At) U*, (21) 

where U = (p, u,p) and the operators Ll(Ai), L2{At), and L3(At) indicate 
the discrete propagators for the corresponding integration steps (the convec- 
tion, the additional, and the acoustic steps, respectively). These equations can 
be summarized as 

U"+^ = L3{At)L2{At)Ll{At) U". (22) 

If we factor L2{At) and L3{At) into some identical pieces, we get 

U"+' = [L3{At/m3)r'^[L2{At/m2)r^Ll{At) U", (23) 



where m2 and m3 are positive integers. This equation means that, after solving 
the convection parts with At, the additional part is solved m2 times with 
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At/m2, and subsequently the acoustic parts are solved m3 times with At/m3. 
If necessary, we divide L2 into two parts as L2 — L2^'^^L2^^\ and modify Eq. 
(23) as 

U"+' = [L3(At/m3) L2^^\At / mS)]"''^ [L2^'\At/m2)r^Ll{At) U" (24) 



or 



U"+' = [L3{At/m3) L2^^\At/m3) [L2^^\At / {m2 m^WY"^ Ll{At) U". 

(25) 

Such a treatment of the additional parts is necessary, e.g., when the surface- 
tension term exists. (The surface-tension part, represented by L2^'^\ should 
be solved together with the acoustic parts because the surface tension always 
needs to balance with the pressure jump at the interface.) Some other forms 
of factorization adjusted for a governing system can be employed. 

We determine the fundamental time interval. At, by 

. ( . h ^ h \ 

At = mm cl c2 — — , 

y max(|u|j max(G5) / 



where the parameters are typically set to cl = 0.2 and c2 = 10. m3 is deter- 
mined so that, e.g., the following condition is satisfied: 



At 

< mm c3- 
^3 \ maxi 



h 



a 



Si+l,j 



CsLj ) 



, c3- 



h 



max 



C. 



At], 



where the spatial difference of sound speed is used as a criterion. If the surface- 
tension term exists, the stability condition needed to solve it [32,33], 



(26) 




should be taken into consideration to determine m3, where pa and p5 are the 
densities of materials on different sides of the interface, and a is the surface- 
tension coefficient. Also, m2 is determined based on the stability condition of 
a method for the viscous term or others. 
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4 Stabilization of the derivative advancement 



At the interfaces where the value or derivatives of physical quantities are in 
general discontinuous, the use of the conventional centered schemes (14)-(16) 
for updating the spatial derivatives is not suitable, providing less accurate 
results, as will be demonstrated below using a numerical example. To modify 
the schemes, we recall here the extrapolation concept investigated in Part 
I [19] and by others [34,35,36,37]. As has been pointed out in Ref. [19], an 
interpolation or differencing across the phase boundary may not physically 
valid, and gives rise to serious numerical errors. In the following, those schemes 
are modified by adapting a simple extrapolation. 



Equation (14) can be rewritten as 

1 / c^j+ij ~ d. 



1,3 



h 



h 



(27) 



1,3 



n+1 



i,J' 



We introduce here a switching parameter H, defined as 



1 for 



otherwise. 



m+l 



>0, 



(28) 



Hij+l/2 



1 for <l>ltl,-<f>lt'>0, 
otherwise. 



(29) 



where is the ID function (a density function or a level set function) defined 
as > in the region occupied by a material and < elsewhere, and is 
updated by solving 

^ + u • V0 = 0. (30) 



Using this, we modify Eq. (27) as 

Q rn+l f) f* _ ^ ( TT dj+ij — djj djj — dj-ij ^ 

OxJij - ()xh,j - 2 1^^+1/2 J + I ■ (31) 

This modification vanishes the derivative ((ij+ij — (ijj)//i when 4>i'+ij-<f>ij^ < 0, 
resulting in the advancement of dxfij using only the values of an identical 
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material. This procedure corresponds to an extrapolation such as dj+ij = dij. 
Equation (15) can be modified in the same manner, as 

Cj/Jij - Cy/ij = 2 1 ^^.^+1/2 ^ I ■ (•^2) 



The scheme (16) for the cross derivative can be rewritten as 



dxyfij^ - dxyfij — ~rr^[{d-i+i,j+i - dij) - {di-ij+i - dij 



We modify this as 



-(di+ij_i - dij) + ((ii_ij_i - dij)]. (33) 



><[-f^i+l/2,i+l/2(c^i+l,j+l — dij) — Hi_i/2,j+l/2{di-l,j+l - dij) 

-Hi+i/2,j-i/2{di+i,j-i - dij) + Hi_i/2,j-i/2{di-ij-i - dij)], (34) 



where 



-f^i±l/2J±l/2 



1 for 0a,±i-C>O, 
otherwise. 



(35) 



The extrapolations used above are merely rough compared with those used 
to solve the convection parts [19,36]; these, however, might be sufficient, be- 
cause the derivatives require lower accuracy than those required for the non- 
derivatives represented here by /. 

Let us perform a numerical test to demonstrate the effectiveness of the above 
modification. A one-dimensional nonlinear sound propagation is solved based 
on Eqs. (4)-(6). The convection and other components are neglected for sim- 
plicity. The initial condition is 



p{x,0) = 




for X < 0.5, 
elsewhere. 



u{x,0)^0, p(x,0) = l, 

X e [0, oo). 
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and the boundary condition is 

p(0,t) = 1.1 - O.lcosa;ii and du{0,t)/dx = Q, 



where the angular frequency cui is determined so that the wavelength of the 
emitted sound wave is 0.225. The initial values of all derivatives are zero. The 
square of sound speed is determined by 

7(p + 3172.04)/pfor x < 0.5, 
lAp/p elsewhere . 




Other parameters are At = 2 x 10~^ and h = 5 x 10~^. (Under this condition, 
K,c is about 0.6 for x > 0.5 and is about 0.15 elsewhere.) The ID function is 
set to be a color function, 

1 for X < 0.5, 
—1 elsewhere. 



The matrix equation given by discretizing Eq. (18) is solved by the Red-Black 
Gauss-Seidel method with the convergence criterion of Eres < 10~^, where 



n+l,{k) _ n+l,(fe-l) 
Pi Pi 



and k denotes the number of the Gauss-Seidel iteration. Figures 1(a) and 1(b) 
show the calculated pressure and density distributions at t = 1.602 x 10~^. 
Also, Figs. 1(c) and 1(d) show the density gradients at the same t, calculated 
using the conventional and the modified schemes. (In the acoustic step, the 
derivatives are not used to update the non-derivatives; thus, the results shown 
in Figs. 1(a) and 1(b) are irrespective of the schemes for the derivatives.) As 
can be clearly seen, the result obtained by the conventional scheme has a 
strong overshoot at the interface, while the result obtained by the modified 
scheme is smooth. The result given using the modified scheme with At /A and 
h/A (the dashed lines) is similar to the latter result. These results can roughly 
prove the validity of the modified scheme, although more detailed discussions 
should be necessary. 



5 Application to the compressible Navier-Stokes equations 



In this section, as a summary of methods, we illustrate how to apply the 
methods discussed in this series to the compressible Navier-Stokes equations 
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with a surface tension term: 



l^ + u- Vp = -pV-u, (36) 



+ u-Vu = ^ + - 2V- (//T) - V(//V • u) + ^, (37) 

or p p \ 6 / P 

|^ + u-Vp=-pC,2V-u, (38) 

where T is the deformation tensor, p is the viscosity coefficient, and F^^ indi- 
cates the surface tension as a volume force. We factor them into the following 
four systems: 

Convection system: 

( dp du dv 



Viscous system: 



Surface-tension system: 

{dp_ ^ 9u^F^ ^ = n 
\at ~ ' dt~ p' dt~ ' 



Acoustic system: 

jdp ^ du Vp dp 



Successively solving these four systems by the MTS technique completes one 
step of the time integration. Among the convection equations, one in terms of 
the density is solved by the hybrid interpolation-extrapolation method intro- 
duced in Part 1, and the others are solved by the conventional CIP (the "Type 
B" multidimensional formula [26] is adapted). The time integration of (f) is 
performed together with the convection system. The viscous system is solved 
by the conventional second-order centered finite differencing, and the surface- 
tension system is solved by the CSF model [32], coupled with the smoothing 
procedure introduced in Appendix B. The acoustic system is solved by the 
modified CUP method reviewed in Sec. 2.2. 
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6 Application results 



This section presents the apphcation results for the multiphase flow problems, 
given using the present methods. Compressible and incompressible Kelvin- 
Helmholtz instabilities and multibubble dynamics in an acoustic field were 
selected as application examples. 

Compressible and incompressible Kelvin-Helmholtz instabilities. It 

is well known that an interface between fluids of different densities is un- 
stable. When the velocity field around the interface is perturbed, the inter- 
face forms a complicated structure [38]. This instability, called the Kelvin- 
Helmholtz instabihty or the Rayleigh- Taylor instability, has been the sub- 
ject of quite a number of theoretical, experimental, and numerical studies 
[38,39,40,41,42,43,44,45], and even the CIP have been employed to investi- 
gate it [41,46]. Using this example, the effectiveness of the extrapolation tech- 
nique for the convection parts and that of the MTS technique are demon- 
strated. The initial arrangement of materials is shown in Fig. 2. The density 
of the heavy material (one for x < 0) is ph — I, and that of the light one 
pi = 1/3, and the initial pressure is po{x. y. t = 0) = 1. The initial velocity 
field [u{x,y,t = 0) = {uo{x,y),Vo{x,y))] is set to 

Uo{x,y) = sgn(a;)Mmaxexp(-27r \x\/\) sm{27iy/X), 

vo{x,y) = tiinaxexp(-27r |x| /A) cos(27r7//A), 

(a;,!/) e [-1.25,1.25] X [0,0.5], 

which satisfies V-u = 0, where Wmax is the maximum velocity and A (= 1) is the 
wavelength of the perturbation. The non-dimensional parameter X-s/phPo/ p is 
fixed to 1400, and p is assumed to be constant. The surface tension and the 
gravity are neglected. The slip boundary condition is adapted to the bound- 
aries of y = and y = 0.5, and the Neumann boundary condition {d/dn = 0) 
to a; = —1.25 and x = 1.25. Figures 3(a)-3(d) show the time sequence of p 
for Umax — 0.7, C| — lAp/p (the Mach number under this condition is about 
0.6), and h (= Ax = Ay) — 0.5/60, at t — 4. Here, the linear monotone 
function [19] is used for the extrapolation, and Eq. (23) is used for the MTS 
integration. The parameters for the time intervals are cl = 0.5, c2 = oo, and 
c3 = 0.1. (The inner iteration number for the acoustic parts, m3, is 7 ~ 13 
in this case.) These results present sharp descriptions of the density interface, 
while the numerical diffusion and oscillation can be observed in the result 
given when the convection parts are solved by the conventional CIP (Figs. 
3(e) and 4). Figures 5(a)-5(d) show = surfaces at t = 4 for different c3. 
(For c3 = oo, the MTS integration is disabled, that is, a single time step is 
used.) As is clearly seen, decreasing c3 changes the result. A similar result to 
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that for c3 = 0.1 can be given, even using a smaller At (see Figs. 5(e) and 
5(f), which show the results for cl = 0.2 and different c3). These results reveal 
that integrating the acoustic parts with a smaller time interval improves the 
accuracy of the solution, even if the time interval for the convection parts 
holds. 

While the above example can be solved by an explicit scheme without any 
difficulty, the next example is problematical. Figure 6 shows the results given 
by artificially resetting the sound speed to C| = 1000 x lAp/p for only the 
heavy material. This setting results in a much greater sound speed than the 
flow velocity, and also results in the coexistence of materials of greatly differ- 
ent compressibilities. The parameters for the time integration are cl = 0.25, 
c2 = oo, and c3 = 5.0. (Under this setting, mS = 3 ~ 5 and ft;c/m3 pa 4 ~ 5.) 
This result is not changed noticeably by decreasing the time interval, proving 
the applicability of the present method to a high kq condition and proving 
its robustness. Using this example, we perform here a convergency test. Fig- 
ures 7(a)-7(c) show the interfaces at t = 4 given using the same parameters 
except for the grid width set to h = ho, ho/ 1.5, or ho/2, respectively, where 
ho = 0.5/60. The refined results are in good agreement with that for h = ho- 
In contrast, the result given by solving the convection parts with only the in- 
terpolation (i.e. by the conventional CIP) changes as the computational grids 
are refined (see Figs. 7(d)-7(f)); the convergency of the solution is obviously 
lower than that obtained by the hybrid method. These results show that the 
improvement given in Part I accelerates the convergency, and that the hybrid 
scheme can accurately describe such deformable interfaces. 

Two-bubble dynamics in an acoustic field. When a sound wave is ap- 
plied, a bubble immersed in a liquid begins volume oscillation. When other 
bubbles exist, they interact acoustically with each other, resulting in the 
changes in the oscillation amplitude and phase and the effective resonance 
frequencies [47,48,49]. Moreover, it is known that an acoustic interaction force 
called the secondary Bjerknes force acts between such pulsating bubbles [50,51]. 
The force is attractive when the bubbles pulsate in-phase, and is repulsive 
otherwise. Such effects resulting from the radiative interaction between bub- 
bles have been important subjects not only for physicists and mechanical and 
acoustical engineers, but also for medical engineers and chemists [52,53,54,55,56], 
and novel insights regarding these effects have been provided even in very re- 
cent years [57,58,59,60,61]. (In Refs. [59,60], for example, the author discovered 
that in multibubble cases, the phase shift of the bubbles' pulsation can take 
place not only at their natural frequencies, but also at some other driving 
frequencies. In a more recent paper [61], it was found that the latter char- 
acteristic frequencies cause the sign reversal of the secondary Bjerknes force 
[51].) 

Solving this problem numerically has some interesting points. (1) The gas in- 
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side the bubbles and the hquid surrounding them have quite different densities 
and compressibihties. (2) The compressibihty of the gas is not neghgible in 
principle, meaning that a solver for incompressible flows is not applicable. (3) 
Especially for small bubbles, the viscosity of the surrounding liquid and the 
surface tension are not neghgible. (4) Many different time scales exist, deter- 
mined by the flow velocity, the sound speed, the surface tension, the viscosity, 
the bubbles' natural frequencies, and the frequency of an external sound. (5) 
The topology of the interfaces changes when the bubbles coalesce as a result 
of their radiative interaction. Based on these interesting factors, we consider 
this problem to be a very good example for demonstrating the abilities of our 
methods. 

The axisymmetric coordinate (r, z) is selected for the computational domain, 
and the centers of the bubbles are located on the central axis. Bubbles 1 and 
2, whose equilibrium radii are Rio — 5 //m and R20 — 9 //m, respectively, are 
filled with a gas of the constant specific heat ratio 7 = 1.33. The equilibrium 
density and the viscosity coefficient are 1.226 kg/m^ and 1.78 x 10^^ Pa s for 
the gas phase, and 1000 kg/m^ and 1.137 x 10~^ Pa s for the liquid phase. 
The surface tension coefficient is cr = 7.28 x 10~^ Pa m, and the equihbrium 
pressure of the hquid is Pq — 101.3 kPa. The initial pressures inside the bubbles 
are thus determined by Pq + 2a/Rjo (for j = 1 and 2), in order to balance 
with the surface tension, where 2a / Riq ^ 0.29Po and 2a/R2o ~ O.I6P0. The 
sound speed is set to C| = jp/p for the gas, and to C| = 7(p + 3172.04Po)/p 
for the liquid. The external sound is applied as the boundary condition to the 
pressure by p(r) = Po(l + 0.3sina;i), where F denotes the boundaries of the 
computational domain except for r = 0, and a; is the angular frequency of the 
sound. The boundary condition for the velocity is free. The initial distance 
between the centers of the bubbles [D{t = 0)] is fixed to 20 /im. The grid 
width assumed to be constant is set to h = (5/20) /xm, and the total number 
of grids is 100 x 310. Pa+Pb in Eq. (26) is set to 1000 kg/m^ (= the equihbrium 
density of the liquid), because the density of the gas is negligible, and that of 
the liquid can be assumed to be almost constant in time. 

Figure 8 shows the bubbles' mean radii and mass centers calculated for u; = Uio 
and cu — a;io/1.8 (~ UJ20), where cujo is the linear, monopole natural frequency 
of bubble j determined by 

u;jo = ^[37Po + (37 - l)2a/Rjo] / pR%, 

and the parameters for the MTS integration were set to cl = 0.2, c2 = 20, 

and c3 = 2, resulting in m3 ~ 8. Also, Figs. 9 and 10 display the interfaces 
at selected times. We can observe the repulsion of the bubbles when u = uio, 
while the attraction when uj = a;io/1.8. These results are validated by the 
theory presented in Ref. [51], and reexamined in Refs. [62,61]. The theory 
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determines the sign of the interaction force (¥23) by 



sgn(F2i3) = sgn(cos(0i - ^2)), 



(39) 



where 



01 



tan-^ (Bi/Ai) e [0, 2n] 



with 
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HiF + M2G 



HiG - M2F 
F^ + G^ 



F 



R10R20 4 



2, 



G 



L1M2 + L2M1, Hi^L2 + 




D 



{ufo-UJ^), L2^{UJ2Q-Uj'^), 



Ml = 6100, M2 = 6200, 

where 6j, set to {Afii/ piR'jq) + {u'^RjQ/ci) in the present study, is the damp- 
ing coefficient [63], /x;, p;, and q are the viscosity coefficient, the equihbrium 
density and the sound speed, respectively, of the hquid, and exchanging 1 and 
2 (or 10 and 20) in the subscripts of these equations yields the expression 
for 02- The positive sign of Eq. (39) indicates attraction, while the negative 
indicates repulsion. Figure 11 shows sgn(F2B) as a function of u ioi D = 20 
/^m, revealing that the bubbles repel each other when uj = cuio, while they at- 
tract when CO — 6i;io/1.8. This theoretical prediction is well reproduced in the 
numerical results. (More detailed discussions from the physical viewpoint are 
given in Ref. [64].) Meanwhile, in the result for cu — cuio/l.S, we can observe 
the coalescence of the bubbles as a result of the attraction. This result proves 
the ability of the present methods to treat the topology change. 

Figure 12 shows At normalized by cl /i/ max(|u|), c2 h/ max{Cs), or Atst, 
as functions of time. This figure reveals that the fundamental time interval 
was mainly determined by the sound speed, whereas the flow velocity was 
responsible when a rapid motion of the interfaces, caused by the coalescence, 
occurred. 

Using the example for uj = Uio, we investigate here the effectiveness of the 
MTS integration on this problem. Figure 13 shows the bubbles' mean radii as 
functions of time, for cl = 0.2 and different c2 and c3. Large differences cannot 
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be seen between the results for c3 < 4, while the results for c3 = 8 (m3 ~ 2) 
are obviously inaccurate. This result proves that the MTS integration has 
obviously contributed to the accurate solutions given for c3 < 4. 

During the computation for oj = cJio/1.8 introduced above, we sometimes 
observed negative p near the interfaces. In such a case, we adopted p = 
max(p, ec-Po), where Ec is set to 10"^; the numerical results were insensitive 
to the magnitude of this parameter. This problem would be overcome by em- 
ploying the oscillation-free variants of the CIP [65,66,67,68]. 



7 Conclusions 

In this series of articles, we have proposed an improved unified solver for com- 
pressible and incompressible fluids involving free surfaces, based on the CIP- 
CUP method, by adapting several improvements and modifications. (The most 
significant one given here in Part II is the adaptation of the MTS integration 
technique, which makes the determination of the time interval very flexible.) 
High accuracy and excellent robustness of the improved methods have been 
demonstrated by using examples of free-surface flows that contain both com- 
pressible and incompressible materials. The present methods, however, face 
the following challenges: 

1. Optimizing the extrapolation function for the convection parts 

In Part I of this series, we proposed flve kinds of extrapolation functions used 
to solve the convection of the density. Although a concrete discussion has not 
been provided in the present paper, we observed that the most accurate result 
for a different problem was achieved by a different extrapolation function. 
The optimization of the function based on some criterion is, therefore, sorely 
expected. 

2. Optimizing the time intervals 

Although in the present study we have empirically determined the parame- 
ters for the MTS integration, an automatic, optimized determination should 
be useful in a practical application. The optimal time interval for the acous- 
tic parts might be determined by a criterion based on the compressibility of 
fluids. If, for example, materials can be considered to be almost completely 
incompressible, as has been well known previously, one can use a time interval 
of an infinite CFL number with respect to the sound speed, whereas when 
a compressible material exists, as has been demonstrated in this paper, one 
needs a time interval depending on the sound speed. This means that the 



18 



maximum values of c2 and c3 sufficient for accurate computations are rightly 
dependent on the compressibility. 

The above subjects will be addressed in a future paper. 
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APPENDIX A: Averaging at phase boundary 



To discretize the acoustic, viscous, and surface tension terms on tiie staggered 
grids, we need the values of p at the velocity positions, i.e., Pi+i/2,j and Pij+i/2- 
In this appendix, we briefly discuss how to estimate them. 

As discussed in Part I, the phase boundary is recognized using the zero level set 
of the ID function (Identification Function; the color or the level set function), 
0, and materials are identified using the sign of the function. If ^j+i j- -^j^j < is 
true, it is recognized that the call between (x^+ij, and {xij, i/ij) contains 

an interface. In such a cell, we estimate the densities by the following VOF 
[69]-like procedure, which implies the weighted linear interpolation: 



Pi+l/2,j 





* 1 

Pi+l,j ' 




Pi,3 
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<0, 
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(41) 



where 0"+^ (= (j)* = 0**. Note that the position of the interface is changed 
only when solving the convection parts.) is the ID function after solving the 
convection part. In the remaining cells, we use the simple average, 

* _ Pi+l,j + Pi,j 
Pi+l/2,j — 2 ' 

* _ Pi,j+l ~^ PiJ 
Pi,j+l/2 - 9 > 



which has been used in the conventional CIP algorithm [18] , or in others. 

The above averaging scheme is also used to estimate the values of some other 
quantities at the velocity positions, such as the viscosity coefficient. 
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APPENDIX B: On the surface tension term 



This appendix concerns how to calculate the surface tension term in our code. 

In the model named the CSF (continuum surface force), proposed by Brackbill 
et al [32] , the surface tension as a volume force is represented by 



is the mean radius of the surface curvature, and 6 denotes the color function 
defined as = 1 for a material, and = otherwise. The singularity appearing 
in V6' at the interface is mollified by some smoothing procedures [32,33,70]. 
In the present study, the smoothed 9 is constructed as follows: 

Step 1: The non-smoothed color function is made from the ID function, which 
is an arbitrary function whose = surface represents the interface [19,42], 
by setting 9ij — 1 for ^j^^ > and 9ij — for (f)ij < 0. (0 always has a 
non-zero value at the gird points [19].) 

Step 2: The values of the color function at the velocity positions, 6'i+i/2j and 
9i,j+i/2, are determined by the averaging procedure described in Appendix A. 

Step 3: Using the 4-point simple average, the weakly smoothed values of the 
color function at the original position is obtained as 

9ij = 0.25(6lj+i/2j + 6'i_i/2,j + 9ij+i/2 + 6'jj_i/2). 
Step 4: Using the following 5-point smoother, 9 is smoothed further: 



where m (= 0, 1, • • • , M — 1) is the number of iteration, M is the total number 
of iteration, and £ is a small positive value. 

Steps 2 and 3 are necessary to approximately take into account the phase 
property of the interface, which is not contained in 9 given at Step 1, but is 



(42) 



where a is the surface tension coefficient, k defined as 




^r'^ = (1 - s)9if + 4(^s, + + 9^, + 9^d, 
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latent between the grid points, e and M are typically set to £ = 0.15 and 
M = 20. 

Replacing 9 with 9'^^\ we calculate the surface-tension term by the second- 
order centered finite differencing on the staggered grids, as is done in Ref. 
[32]. 
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Fig. 1. Nonlinear propagation of sound wave. Pressure (a), density (b), density 
gradient calculated by the present method (c), and by the conventional method (d), 
at t = 1.602 X 10~^. The dashed lines denote the result for At/4 and h/A. 
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Fig. 2. Initial material arrangement for Kelvin-Helmholtz instabilities. 
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Fig. 3. Results for compressible Kelvin-Helmholtz instability by the hybrid method 
(a-d) and by the conventional CIP (e). 
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Fig. 4. Density profiles on a; = 1.5 at t = 4, by the hybrid method (left) and by the 
conventional CIP (right). 
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Fig. 6. Kelvin-Helmholtz instability in a composite flow of compressible and incom- 
pressible fluids. 
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Fig. 7. Convergency test. Results at t = 4 given by the hybrid method (a-c) and 
by the conventional CIP (d-f), for h = ho, ho/1.5, and ho/2. 
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Fig. 8. Bubbles' mean radii and positions (the lower one is for the smaller bubble) 
for u> = uio and u = wio/1.8 as functions of time. The coalescence of the bubbles 
occurred at the time when the number of the lines became one. 
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Fig. 11. Theoretical result for the sign of the secondary Bjerknes force for D = 20 
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Fig. 12. Time intervals normalized by cl /i/ max(|u|) {Kl), c2h/ msix(Cs) (^2), 
or Atst {K3) as functions of time. 
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Fig. 13. Bubbles' mean radii for u = uio as functions of time, for cl 
different c2 and c3. 
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